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1.  Introduction 


Parabolic  approximations  to  the  reduced  wave  equation  have  been  used 
extensively  in  acoustic  propagation  since  the  early  1970’s  [1]  and  in 
electromagnetics  even  earlier  [2].  The  advantage  of  a  parabolic  equation  (PE) 
is  that  the  solution  can  be  marched  forward  in  range;  the  field  at  a  given  range 
Vq  is  not  dependant  on  the  field  at  ranges  r  >  Vq.  In  contrast,  the  reduced  wave 
equation  is  elliptic,  so  the  field  at  range  Kq  is  dependant  on  the  field  at  all  other 
ranges.  This  requires  solving  a  large  set  of  simultaneous  equations,  a  much 
more  demanding  problem  computationally. 

Various  PEs  have  been  implemented  numerically  using  finite-differences, 
finite-elements,  and  a  Fast  Fourier  Transform  (FFT)  based  method  known  as  the 
split-step  PE.  The  split-step  PE  has  a  significant  advantage  in  that  the  range 
step  typically  is  on  the  order  of  tens  of  wavelengths  as  opposed  to  tenths  of  a 
wavelength  for  the  finite-element  and  finite-difference  models.  This  makes  the 
split-step  much  faster,  usually  by  at  least  one  order  of  magnitude. 

The  underwater  acoustic  community  was  the  first  to  apply  the  PE,  including  the 
split-step,  to  the  acoustic  propagation  problem.  The  biggest  problem  faced  in 
adapting  the  underwater  models  to  the  problem  of  atmospheric  propagation  is 
accommodation  of  the  air-ground  boundary.  The  Green’s  function  parabolic 
equation  (GFPE)  [3]  is  the  first  split-step  PE  to  include  a  complex-impedance 
surface.  The  purpose  of  this  report  is  to  present  a  theoretical  development  of 
PEs  in  general,  culminating  in  a  detailed  description  of  the  GFPE.  It  brings 
together  theory  from  a  number  of  sources  in  a  cohesive  manner  with  consistent 
notation.  It  is  the  first  in  a  series  of  three  reports  on  the  GFPE.  The  second 
report,  A  Sensitivity  Study  of  the  Green ’s  Function  Parabolic  Equation, 
examines  the  sensitivity  of  the  GFPE  to  the  input  parameters  and  presents  an 
approach  to  the  automated  selection  of  several  of  the  parameters  based  on  the 
sound  speed  profile.  The  third  report,  A  Users  Guide  to  the  Green ’s  Function 
Parabolic  Equation,  gives  a  brief  review  of  the  first  two  and  describes  the  user 
interface  from  a  more  user  oriented,  less  theoretical  point  of  view. 

Section  2  derives  the  reduced  wave  equation,  the  starting  point  for  parabolic 
approximations.  Section  3  presents  the  original  narrow-angle  PE,  followed  by 
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a  discussion  of  the  general  split-step  approach  and  some  wide-angle 
approximations.  Section  4  applies  the  split- step  technique  to  the  problem  of  the 
air-ground  interface,  leading  to  the  GFPE.  Section  5  gives  examples  of  the 
output  of  the  GFPE,  along  with  some  comparisons  to  another  outdoor  sound 
propagation  model,  the  Fast  Field  Program  (FFP). 


2.  The  Reduced  Wave  Equation  for  Acoustic  Propagation 


Begin  with  two  conservation  equations  from  fluid  mechanics:  the  conservation 
of  mass, 

^  +  V-(pV)  =  0,  (1) 

dt 

and  the  conservation  of  momentum,  neglecting  viscosity*, 

— (pV)  +  VP  =  0  (2) 

dt 


where 

p  =  the  density  of  the  fluid 
P  =  the  pressure 
V  =  the  fluid  velocity. 

Each  is,  in  general,  a  function  of  position  and  time.  Furthermore,  P  is  assumed 
to  have  some  functional  dependance  on  density  and  entropy  S, 

P  =  />(p,«.  (3) 

Acoustic  disturbances  modeled  by  PEs  involve  very  small  fluctuations  from 
ambient  levels,  so  the  following  assumptions  are  made: 

P  =  Po  +  p',  P  =  Pq  +  P', 


*Eor  the  acoustic  phenomenon  discussed  here,  the  effects  of  viscosity  are  negligible. 
In  a  more  general  setting,  viscosity  may  be  important. 
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with 


(5) 

Po  Po 


The  primed  variables  are  the  acoustic  fluctuations,  while  the  subscripts  indicate 
ambient  quantities.  Note  that  the  total  velocity  is  assumed  to  be  small.  The 
PEs  discussed  here  assume  an  ambient  wind  velocity  of  zero.  Inserting 
equation  (4)  into  the  conservation  laws  and  neglecting  products  of  small 
quantities  yields  the  linearized  conservation  laws 

+  p„V-P  =  0  (6) 


and 


p„f  .  V/  =  0. 


(7) 


Small  amplitude  sound  propagation  is  generally  isentropic,  so  that 

P  =  P(p).  (8) 


This  can  be  linearized  around  the  point  pg, 


c^p'. 


c"  = 


ap(p)  I 

3p  p  ■  Po 


(9) 


Equations  (6),  (7),  and  (9)  can  then  be  solved  for  p '  resulting  in  the  acoustic 
wave  equation 

=  0,  (10) 


where 

c  =  c(j^  =  the  sound  speed. 
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The  next  step  is  to  assume  a  source  periodic  in  time,  with  angular  frequency  co 
This  can  be  expected  to  give  pressure  fluctuations  of  the  form 

so  that  p  satisfies  the  reduced  wave  equation 

VV  +  k^p  =  0,  (12) 


where 

k  —  co/c. 

Using  cylindrical  coordinates  (figure  1)  and  assuming  the  pressure  field  is 
independent  of  0  leads  to 

^  ^  =  0.  (13) 

dr^  r  dr  dz^ 
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It  is  well  known  that  a  cylindrical  wave  decreases  in  amplitude  as  To 
separate  this  cylindrical  spreading  loss  from  other  effects,  define  u  by 


p{r,z)  -  r  ^u{r,z). 
Then  u  satisfies  the  reduced  wave  equation 


dr^ 


/j2  In  /n 

+  +  )u  =  0. 

dz^  4r^ 


Finally,  assume  (krf>  1,  which  gives  rise  to 


dr^ 


(14) 

(15) 

(16) 


Equation  (16)  is  the  reduced  wave  equation  describing  acoustic  propagation  far 
from  the  source.  It  is  the  starting  point  for  the  parabolic  approximations  of  the 
following  sections. 


3.  An  Overview  of  Parabolic  Approximations 


In  the  development  that  follows,  the  wavenumber  k  is  assumed  independent  of 
r,  but  is  allowed  to  vary  with  z.  In  practice,  parabolic  models  are  used  in 
weakly  range  dependant  environments.  In  that  case,  some  of  the  results  below 
become  approximations  rather  than  exact  representations.  This  is  pointed  out 
where  it  occurs.  It  turns  out  that  the  approximations  are  generally  good,  and 
the  models  yield  good  results  if  the  range  dependance  is  not  too  drastic. 

3.1  The  Narrow- Angle  Approximation 

The  first  use  of  a  PE  for  acoustic  propagation  was  that  of 
Tappert  and  Hardin  [1],  who  invoked  the  restriction  to  narrow  propagation 
angles.  Their  model  is  presented  in  this  section. 

Begin  by  expressing  u  as  the  product  of  a  horizontal  plane  wave  of 
wavenumber  kQ  and  an  envelope  'P, 

u{r,z)  = 


where 


k(,  —  &Icq  for  some  reference  sound  speed  Cq. 


Note  that  'P  =  |  w  |  only  when  k  =  kg,  which  can  only  occur  in  a  homogenous 
medium.  Substitution  of  this  into  the  reduced  wave  equation  (16)  yields 


dr^ 


+  lik^ 


ay 

dr 


a^y 

az" 


+  ik^  -  -  0. 


(18) 


This  can  easily  be  made  into  a  PE  by  omitting  the  second  partial  derivative  of 
y  with  respect  to  r, 

lik^—  +  ^  -  kl)^  =  0.  (19) 

"  ar 
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But  when,  if  ever,  is  this  justified?  A  sufficient  condition  is 


dr^ 


< 


\k. 


ar 

dr 


(20) 


For  a  homogenous  atmosphere  of  wavenumber  kg,  this  can  be  interpreted  via 
the  Fourier  Transform  pair  at  a  given  r  = 


iLz  77 

e  ^  dk„ 


(21) 


and 


oo 

uirjc^  =  j  e  '’^^'^u(r,z)dz. 

—  oo 


(22) 


Equation  (21)  can  be  viewed  as  a  superposition  of  plane  waves  of 
wavenumber  kg. 


oo 

u(rQ+Ar,z)  =  dk^. 


(23) 


where 


-  ko~ 

Ar  =  0. 


For  nonzero  Ar,  this  represents  propagation  of  to  r  =  +  Ar  in  a 

homogenous  medium  of  wavenumber  kg.  u(rQ,k^)  is  the  amplitude  of  the  plane 
wave  component  whose  angle  of  propagation  0  satisfies 


cos  6  = 


(24) 
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See  figure  2.  In  terms  of  the  envelope,  this  becomes 


¥(ro  +  Ar,z)  = 


oo 

—  oo 


(25) 


propagation  of  plane  wave 
components. 


Applying  the  sufficient  condition  of  equation  (20)  with  respect  to  Ar  gives  rise 
to  \k^-  ko  \  <ko,  or 

COS0  «  1,  (26) 

hence,  the  name  narrow-angle  approximation.  In  the  case  of  a  nonhomogenous 
atmosphere,  the  wavenumber  variations  will  introduce  wavefront  perturbations 
so  equation  (25)  does  not  apply.  However,  in  most  cases,  the  perturbations  are 
very  small  and  the  wavefronts  are  still  locally  approximately  planar.  Thus, 
equation  (20)  is  still  roughly  satisfied  by  the  narrow-angle  condition  of 
equation  (26).  More  is  said  about  this  in  sections  3.3  and  3.4. 

3.2  Operators  and  the  Split-Step  PE 

A  more  general  approach  to  the  PE  involves  approximations  of  a  certain 
differential  operator,  which  will  be  considered  next.  But  first,  it  is  interesting 
to  observe  that  the  term  PE  is  really  a  misnomer,  because  most  of  the  operator 
approximations  result  in  equations  that  are  not  parabolic.  However,  this  paper 
follows  custom  and  continues  to  use  the  name  parabolic  for  equations  in  this 
class.  Also,  be  aware  that  frequent  use  is  made  of  various  representations  for 
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functions  of  an  operator,  without  mathematical  justification.  The  formal 
manipulations  are  well  defined,  and  the  results  are  valid  only  for  operators 
satisfying  certain  unstated  conditions.  No  attempt  is  made  to  prove  the 
operators  involved  satisfy  necessary  or  sufficient  conditions.  The  ultimate 
justification  is  empirical;  the  final  results  agree  well  with  experimental  data  and 
other  models. 

Begin  by  expressing  equation  (16)  in  operator  notation, 

[■4  ^  klQ]u  =  0,  (27) 


where 


Q  =  + 


kl 


(28) 


n  =  k/kg  =  the  acoustic  index  of  refraction. 

Because  k,  and  therefore  n,  is  independent  of  r,  equation  (27)  can  be  factored 
as  follows: 

-  iK-jQ)  U  =  0.  (29) 

or  dr 

The  left  and  right  factors  represent  inwardly  and  outwardly  propagating  waves, 
respectively.  If  k  depends  on  r,  the  partial  derivative  with  respect  to  r  does  not 
commute  with  so  the  factorization  is  only  an  approximation. 

Interest  is  in  the  waves  propagating  outward  from  the  source,  so  only  the  right 
term  is  considered. 


—  =  ikfjQu. 
dr 


(30) 
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The  formal  solution  to  this  equation  is  [4] 


u(r+Ar,z)  =  e 


i^rk^y/Q 


uir,z). 


Note  the  similarity  to  the  ordinary  differential  equation 

dx 

—  -  ax  , 
dt 


which  has  the  solution 


x{t+t^  =  X  (t^. 


The  difference  of  course  is  that  the  exponential  on  the  right  side  of 
equation  (31)  is  a  differential  operator.  The  essence  of  the  split-step  technique 
involves  approximations  of  the  square  root 

!&  se-  (34) 


=  (!+€  +  [ly. 


where 


G  =  n^-1 


kl 


The  general  technique  is  to  linearly  separate,  or  split  the  operators  ^  and  s  as 

y/Q  ^  ^  =  A(\i)  +  B(e), 
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then  write 


^iArk^^/Q  ^  (38) 

The  symmetric  splitting  of  the  exponential  is  exact  only  when  A  and  B 
commute,  which  will  not  occur  if  k  varies  with  z.  However,  it  may  still  prove 
to  be  a  good  approximation,  as  is  discussed  in  the  next  section. 

This  split-step  operator  has  a  simple  interpretation.  A  is  independent  of  s  and 
therefore  n,  so  the  exponential  terms  in  A  represent  propagation  through  a 
homogenous  atmosphere.  All  of  the  variation  in  n  occurs  in  the  exponential 
term  containing  B,  and  5  is  a  multiplicative  operator.  Thus,  the  operators  in 
equation  (38)  represent  propagation  through  a  homogenous  atmosphere  a 
distance  of  Ar/2,  followed  by  a  multiplicative  phase  correction  caused  by 
variations  in  k,  followed  by  another  propagation  through  a  homogenous 
atmosphere  a  distance  of  Ar/2. 

In  practice,  split-step  PEs  are  usually  implemented  using  a  nonsymmetric 
splitting, 

^iArk^s/Q  _  ^iArk^^iArk^B  (39) 

The  field  is  initially  propagated  through  free  space  a  distance  of  Ar/2  and  then 
alternately  phase  corrected  and  stepped  forward  by  Ar.  It  turns  out  that  the 
error  caused  by  noncommutativity  of  A  and  B  is  reduced  by  using  a  symmetric 
splitting  [5],  but  past  the  first  Ar/2  step,  the  symmetric  and  nonsymmetric 

splittings  yield  the  same  result,  because  =e  2  2  “  The  real  issue 

is  the  approximation  of  the  operator  \[Q  . 

3.3  Some  Wide-Angle  Approximations 

Several  approximations  have  been  discussed  in  the  literature  [6,7,8].  One 
consists  of  taking  the  linear  terms  in  the  formal  series  for  the  square  root  of 
1  +x. 
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(40) 


V5  =  =  1  .  |e  . 

Another  approximation,  originally  suggested  by  Tappert  [6],  is 

v/5  =  =  (1  +  +  ie.  (41) 

This  will  be  referred  to  as  the  GFPE  approximation  because  it  is  used  by 
Gilbert  in  the  GFPE  discussed  in  the  next  section. 

The  wide-angle  approximation  of  Feit  and  Fleck  [7]  is 

TG  -  -  (1  +  +  (1  +  G)^-l. 

These  approximations  require  different  restrictions  on  s  and  p,.  Intuitively,  it 
is  necessary  to  consider  the  "magnitude"  of  these  operators.  However,  the  usual 
operator  norm, 

\\A\\  =  sup  \Au\,  (43) 

l«l  =  1 

will  not  work.  Because  p  is  unbounded,  the  supremum  is  oo.  Instead,  as  an 
angle-dependant  measure  of  the  magnitude  of  A,  define 

IHIe  =  \Au\,  (44) 


where  u  is  restricted  to  unit  amplitude  plane  waves  u  -  of  wavenumber 

A:/  =  A:/  +  and  propagation  angle  0  (tan  0  =  kjk^. 

Note  that 

II  pig  =  sin^e,  (45) 

so  p  is  small  exactly  when  0  is.  Also,  |  s  |  ^  1  is  equivalent  to  «  1,  or  A: «  kg. 
One  would  expect  that  the  first  approximation  (equation  (40))  applies  for  small 
p  and  s,  and  these  are  precisely  the  conditions  used  for  the  narrow-angle  PE  of 
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section  3.1.  It  is  not  difficult  to  show  that  equation  (40)  results  in  the 
narrow-angle  PE. 

More  careful  analyses  of  these  approximations  can  be  made  following 
Thompson  and  Chapman  [8].  Define  the  error  operators 

E^=  -  Q  =  +  eq  +  jre  +  ix\  (46) 


ll£^|le  <  ^|6n|(2  +  |6n|)[|6«|(2  +  |6ni)  +  4|cos0-l|],  (50) 

and 

llfp^llg  <  2|6n|  |cos6-l |.  (51) 

Figures  3  through  6  show  the  error  bounds  as  a  function  of  0  and  5n, 
normalized  by 

WQWq  =  (1  +  6nf  -  sm20.  •  (52) 

Figures  3  and  4  vary  0  for  constant  5n  of  0.01  and  0.1,  respectively.  Note  that 
for  0  =  0,  the  wide-angle  approximation  is  exact,  while  the  other  bounds  are 
nonzero.  As  0  increases  beyond  about  20°,  the  error  bound  of  the  narrow-angle 
approximation  begins  to  climb  rapidly.  The  GFPE  and  wide-angle 
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approximations  are  seen  to  be  valid  for  larger  0,  the  difference  between  the  two 
being  the  constant  bias.  Thus,  these  approximations  are  considered  wide  angle. 


Upper  Bound  on  Relative  Error,  5n=0.01 


«  0.04 


(U 

0.02  h 


'  10  20  30  40 

Angle  (degrees) 


Figure  3.  Relative  errors  as  a  function  of  d  for 
6n  =  .01.  Eyv  is  wide  angle,  is  GFPE,  E^  is  narrow 


Upper  Bound  on  Relative  Error,  5n=0.1 
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Figure  4.  Relative  errors  as  a  function  of  9  for  on  = 
Ew  is  wide  angle,  Eg  is  GFPE,  E^  is  narrow  angle. 


Figures  5  and  6  show  the  normalized  errors  as  a  function  of  5n,  for  constant  0 
of  20°  and  40°,  respectively.  The  narrow-angle  approximation  is  seen  to  have 
a  nonzero  bound  for  5n  =  0,  while  the  other  two  are  exact.  5n  =  0.1  represents 
a  sound  speed  variation  of  approximately  30  m/s,  which  is  larger  than  almost 
any  variation  encountered  in  the  surface  boundary  layer.  Thus,  the  wide-angle 
and  GFPE  approximations  are  again  seen  to  be  good  up  to  a  propagation  angle 
of  40°. 


Upper  Bound  on  Relative  Error,  0=20  degrees 


Figure  5.  Relative  errors  as  a  function  of  6n  for  0  =  20. 
is  wide  angle,  is  GFPE,  E^  is  narrow  angle. 


Figure  6.  Relative  errors  as  a  function  of  6n  for  6  =  40. 
is  wide  angle,  Eg  is  GFPE,  E^  is  narrow  angle. 


McDaniel  [9]  studied  the  error  associated  with  the  narrow-angle  approximation 
by  applying  separation  of  variables  to  the  reduced  wave  equation  and  the 
narrow-angle  PE,  equation  (19).  The  normal  modes  for  the  reduced  wave 
equation  in  the  far  field  are 


p„(r  +  Ar,z)  =  ( 


r  +  Ar 


ik„Ar  .  . 


(53) 


while  the  far  field  normal  modes  for  the  PE  are 


p„(r  +  Ar,z)  -  ( 


K  *  K 


r  +  Ar 


Ye 


2Ao 


-Ar 


P„(r,z). 


(54) 


In  both  cases,  k„  is  the  separation  constant,  determined  by  the  boundary 
conditions  on  the  z-dependant  equation.  This  equation  is  identical  for  the 
reduced  wave  equation  and  the  PE.  Thus,  the  PE  is  seen  to  posses  a  phase 
error,  resulting  from  the  wavenumber 
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(55) 


7  2  il 

K  ^  K 

2k,  ’ 


where  k,,  is  the  correct  wavenumber  obtained  from  the  reduced  wave  equation 
solution.  This  results  in  a  phase  velocity  error  in  the  time-dependant  solution 


^  2 
^0  c„ 


1  ^-l 


(56) 


where  =  the  correct  modal  phase  velocity  for  the  wave  equation.  This  error 
is  the  result  of  the  narrow-angle  approximation;  it  is  not  a  numerical  error 
associated  with  any  particular  numerical  implementation  technique.  The 
wide-angle  and  GFPE  approximations  can  be  expected  to  introduce  some  phase 
errors  as  well,  although  quantitative  analyses  of  these  errors  are  more  difficult. 


3.4  Implementation  Via  the  FFT 


Having  considered  the  error  associated  with  the  wide-angle  approximation 
it  remains  to  implement  it  in  a  computationally  fast  algorithm.  The  key  is  to 
observe  that  for  s  =  0,  gives  an  exact  solution  to  the  Helmholtz  equation 
with  constant  wavenumber  kg.  Recall  that  with  the  nonsymmetric  splitting, 

u(r  +  6.r,z)  =  ^ 


where 

A  =  \J\  +  \\.  for 
B  =  \/l  +  e  -  1  =«  -  1  for 

Thus,  when  s  =  0,  (n  =  1) 

u{r  +  Ar,z)  = 

describes  one-way  propagation  through  homogenous  space.  Section  3.1  shows 
that  this  call  be  represented  by 
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00 

M(ro  +  Ar,z)  =  ^  /  ^ 


(59) 


where 


u 


the  Fourier  Transform  of  u  with  respect  to  vertical 
wavenumber  k. 


+  k!  = 


Kg  . 


But  this  also  represents  the  inverse  Fourier  Transform  of  the  quantity 


expressed  as 


ue  V  -  -  _  Thus,  for  the  wide-angle  PE,  the  homogenous  propagation  is 


(60) 


where 

^  =  the  Fourier  Transform  operator  with  respect  to  z 

=  the  associated  inverse  transform. 

Inclusion  of  the  phase  correction  term  for  a  nonhomogenous  atmosphere  yields 

...  iArkf^B  iArkf^A  .  . 

M(r+Ar,z)  ~  e  “  e  “  u(r,z) 

r~2  2  / 

=  /M(r,zO]  • 

This  is  the  wide-angle  split-step  PE.  The  Fourier  Transforms  are  implemented 
via  an  FFT,  which  can  be  computed  with  speed  and  efficiency  using  one  of  the 
widely  available  power-of-two  algorithms.  In  practice,  the  range  steps  are 
several  wavelengths  long,  in  contrast  to  finite-element  and  finite-difference 
techniques  that  use  range  steps  much  less  than  a  wavelength.  The  combination 
of  long  range  steps  and  FFT  is  what  makes  the  split-step  PE  so  fast. 
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Notice  that  no  explicit  mention  of  boundaries  such  as  the  ground  or  water 
surface  has  been  made.  In  general,  integral  transforms  are  applied  so  as  to 
explicitly  represent  the  boundary  conditions.  For  the  present  problem,  that 
would  require  a  transform  with  respect  to  range  r,  so  that  the  boundary 
condition  would  be  translated  to  u(k„z).  Instead,  the  split-step  PE  applies  the 
transform  with  respect  to  z,  so  that  the  boundary  no  longer  exists  in  the 
wavenumber  domain.  Thus,  the  boundary  conditions  can  not  be  accounted  for 
explicitly;  they  actually  become  conditions  on  the  symmetry  of  ^(r,!^). 

A  condition  which  can  be  incorporated  is  that  of  a  pressure  release  surface, 
u(r,0)  =  0.  In  this  case,  the  function  u  can  be  defined  for  z  <  0  by  u(r,z)  =  - 
u(r,-z),  which  forces  the  pressure  release  condition.  Then  the  transform  is 
defined  for  all  z,  eliminating  an  explicit  boundary  condition.  This  is  the 
method  of  images,  modeling  the  reflected  waves  as  direct  waves  propagating 
upward  across  the  z  =  0  boundary.  The  reflection  coefficient  is  R  =-  1. 

Of  course,  the  FFT  imposes  two  boundary  conditions,  because  the  transform 
must  be  finite  in  length.  The  lower  boundary  can  be  treated  as  a  pressure 
release  surface  by  exploiting  the  relationship  between  symmetries  in  the  two 
domains;  the  reflected  spectrum  is  -1  times  the  mirror  image  of  the  direct 
spectrum.  Because  of  the  periodicity  induced  by  the  FFT,  this  will  result  in  a 
pressure  release  surface  at  the  other  boundary  as  well.  The  unwanted 
reflections  from  this  surface  can  be  attenuated  by  introducing  an  artificial 
absorption  layer  near  the  upper  boundary.  Care  must  be  taken  not  to  attenuate 
refracted  sound  which  may  reach  the  detector. 

This  model  works  well  for  underwater  acoustics,  where  it  has  been  used  since 
the  early  1970’s.  The  air-water  interface  is  a  pressure-release  surface  (at  least 
for  smooth  surfaces),  and  the  ocean  bottom  is  generally  a  thick  layer  of 
sediment  that  strongly  attenuates  sound  with  little  reflection.  However,  this  is 
not  the  case  for  outdoor  sound  propagation,  as  the  ground  generally  is  neither 
a  pressure-release  surface  nor  a  thick  attenuating  layer.  Thus,  the  split-step  PE 
needs  a  formulation  that  explicitly  accommodates  more  general  boundary 
conditions. 


4.  Complex  Ground  Impedance  and  the  GFPE 


Reflections  from  the  ground  can  have  a  significant  effect  on  the  propagating 
wavefield.  Multiple  reflections  in  a  downward  refracting  atmosphere  can  result 
in  propagation  to  very  great  distances.  Coherence  effects  between  direct  and 
reflected  wavefields  can  produce  interference  patterns.  To  incorporate  these 
effects  into  a  model,  it  is  necessary  to  accurately  represent  the  ground 
interaction. 

In  many  cases,  the  ground  can  be  modeled  as  a  locally  reacting  surface,  using 
a  complex  surface  impedance  [10].  The  impedance  can  easily  be  incorporated 
into  finite  element  and  finite  difference  PEs,  but  the  wide-angle  split-step  PE 
incorporates  only  a  pressure-release  surface,  as  discussed  above.  The  advantage 
of  the  GFPE,  as  shown  below,  is  that  it  incorporates  a  complex-impedance 
surface,  while  still  allowing  implementation  via  an  EFT  for  speed.  It  is  the 
combination  of  these  two  features  that  distinguishes  the  GFPE  from  previous 
implementations . 

4.1  The  Split-Step  Approximation 

The  starting  point  is  again  the  outgoing  component  of  the  reduced  wave 
equation  (16), 

^  =  i/Lv/Ou,  (62) 

dr 


The  vertical  dependance  on  height  is  characterized  by 


25 


k\z)  =  V  + 

where  the  perturbation  b]c(z)  is  assumed  small, 


(64) 


,2 


1. 


(65) 


In  terms  of  n  and  s,  this  becomes 


1  +  ^ 


=  1  +  G. 


(66) 


The  GFPE  is  a  split-step  PE  with  the  approximation  Oq  discussed  above. 


v/5  =  (1  +  li)"  t  ^  =  (1  +  4"^)”  ^ 

2  kl  dz^  2kl 


(67) 


A  nonsymmetric  splitting  is  used,  giving  rise  to  the  formal  solution 


u(r  *  Ar,z)  = 


ils.rbk^ 

J 


Ar^ox/O 


u{r,z). 


(68) 


where 


A  =  (1 


1 


)"  =  M 


(69) 


As  with  the  wide-angle  split-step  PE,  this  solution  first  propagates  a  distance 
Ar  through  a  homogenous  atmosphere  of  wavenumber  k„,  and  then  applies  a 
phase  correction  for  the  nonhomogeneous  wavenumber  variations.  Note  that  the 
wide-angle  PE  and  GFPE  use  the  same  operator  for  the  homogenous 
propagation.  The  distinguishing  feature  of  the  GFPE  is  its  incorporation  of  a 
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complex  impedance  boundary  condition.  To  do  this,  a  spectral  representation 
is  employed  to  express  as  the  integral  of  a  Green’s  function. 

4.2  The  Spectral  Representation  of  an  Operator 

The  spectral  form  for  a  function  of  an  operator  A  is  [4] 

m  =  fix)  [xl  -  AY'  dx,  (70) 

2%i  c 

where 

C  =  the  boundary  of  an  open  set  containing  the  spectrum  of  A. 

The  spectrum  is  defined  so  that  the  operator  [xl  -  A]  is  invertible  and  bounded 
everywhere  on  C  and  in  the  region  bounded  by  C.  In  the  present  case,  the 

function  is  giving  rise  to 

=  J-<f  -  QY^dx.  (71) 

2ni{, 

The  change  of  variables  x  =  leads  to 

=  if  [k^I  -  liQ\-'Kdk^  (72) 

K' 

This  is  an  operator,  and  its  action  on  u,  may  be  found  by  bringing  u 

inside  the  integral, 

u(r,z)  ^  -  klQy^u(r,z)k^dk^.  (73) 

d 

Observe  that  the  integrand  is  no  longer  operator  valued.  Defining  the 
function  (])  as 

^{K.r,z)  =  [kll  -  klQY^uir,z) 
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leads  to 

u(r^)  =  —Se^^''<i>(k,rX)kdk  (75) 

KiJ 

In  terms  of  the  definition  of  Q,  the  integrand  <{)  satisfies  the  equation 


(-^  +  (76) 

dz^ 

with  =  k^  -  k^.  This  can  be  solved  via  a  Green’s  function  G(k^,z,z'), 

00 

=  f  G(k^^^^)u(r^^)dz', 

0 


where  G  satisfies 

(~  +  =  -6(j  -  z').  (78) 

dz 

All  together, 

00 

u(r;z)  =  G{k^^^')u{r;z')  dz'k^dk^.  (^9) 

0 

The  next  step  is  to  solve  for  the  Green’s  function,  incorporating  the  complex 
impedance  boundary  condition. 

4.3  Incorporation  of  Complex  Ground  Impedance 

The  particular  form  of  G  depends  on  the  choice  of  boundary  conditions  for 
equation  (78).  In  general,  solutions  contain  direct  and  reflected  terms  like 

ik 

e  ^  ,  where  Z{z,z  0  is  a  piecewise  linear  function  such  as  z  +  z '  or  |  z  -  z '  | . 

Thus,  equation  (79)  constructs  the  field  as  a  superposition  of  plane  wave 
components 
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(80) 


^i[k^Ar  *  k^(.z^)] 

The  correct  boundary  conditions  for  equation  (78)  are  a  radiation  condition  for 
large  z  and  the  surface  boundary  condition  at  z  =  0. 

In  general,  the  surface  boundary  condition  can  be  quite  complicated,  describing 
nonlocal  interaction  between  the  air  acoustic  waves  and  the  ground.  A  plane 
wave  impinging  upon  the  surface  at  one  point  can  induce  a  surface  wave  that 
may  interact  with  waves  at  another  location.  In  many  cases,  the  surface  can  be 
modeled  as  locally  reacting.  In  this  case,  plane  waves  incident  upon  the  surface 
are  strongly  refracted,  so  the  transmitted  energy  propagates  nearly  perpendicular 
to  the  surface,  producing  only  local  interaction.  The  boundary  condition  for 
this  type  of  surface  can  be  described  using  surface  impedance. 

The  surface  impedance  is  defined  as 

Z  =  (81) 

v’ 


where 

p’  =  the  acoustic  pressure  at  the  surface 

V„  =  the  outward  component  of  acoustic  velocity  normal  to  the  surface. 

Because  the  pressure  and  velocity  may  not  be  in  phase,  Z  is  generally  complex. 
When  the  surface  is  the  ground,  it  is  referred  to  as  the  complex  ground 
impedance.  There  are  many  models  available  to  estimate  the  ground 
impedance  [10,1 1].  While  many  models  are  based  on  physical  principles,  most 
are  somewhat  phenomenological  or  empirical. 

Boundary  conditions  for  a  differential  equation  such  as  equation  (78)  are 
generally  specified  as 

a +  fcG(r.z)  =  01,  (82) 

OZ 
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where  >  0.  Using  the  linearized  conservation  of  momentum 

equation  (7)  and  the  definitions  of  Z  in  equation  (81),  p  in  equation  (11),  and 
=  co/cq,  the  vertical  momentum  equation  becomes 


dz 


dt 


Po  dp' 


"Po./  _  Po^. 


^  *7  *  U 

-P  -  -iK-ttP- 


(83) 


This  applies  to  the  Green’s  function  G,  which  represents  acoustic  pressure.  In 

Z 

terms  of  the  dimensionless  quantity  Z  =  - ,  the  boundary  condition  is 

Po^ 


^  .  i^G  -  0 

dz  Z„ 


Iz  =  0- 


(84) 


The  resulting  Green’s  function  is 


2k^ 


*  R(.K)e 


ik(z  +  z^y 


], 


(85) 


where 


Rik;)  - 


K^g  +  K 


(86) 


is  the  plane  wave  reflection  coefficient.  Substituting  this  Green’s  function  into 
equation  (79)  gives 


u{r,z)  =  f 

TZl''  '' 


ikjz  -  z'\ 


2k^ 

''\u{r,z')dz'k^dk^. 


(87) 


Implementation  Via  the  FFT 


The  final  goal  is  to  implement  the  GFPE  via  an  FFT,  as  with  the  wide-angle 
split-step  PE  of  section  3.  Thus,  it  is  necessary  to  manipulate  equation  (87) 
into  the  form  of  a  Fourier  Transform  with  respect  to  z,  followed  by  a 


multiplication,  followed  by  an  inverse  transform  with  respect  to  a  vertical 
wavenumber.  The  Green’s  function  is  already  expressed  in  terms  of  z  and 
but  the  contour  integral  is  with  respect  to  horizontal  wavenumber  k^.  Thus,  it 
is  necessary  to  introduce  an  integration  with  respect  to  vertical  wavenumber  and 
perform  the  horizontal  wavenumber  integration. 


Therefore,  consider  the  integral 


J_ 

2k 


/ 


giK(z  -  z') 

-  k: 


Jk, 


(88) 


which  can  be  evaluated  using  the  method  of  residues.  The  standard  technique 
is  to  integrate  around  a  closed  contour  Q  in  the  complex  plane  consisting  of  the 
interval  {-a,d\  along  the  real  axis  and  a  semicircle  of  radius  a  in  the  upper  or 
lower  half  plane  (a  >  0).  The  choice  of  half  plane  is  made  to  ensure  that  the 
integral  along  the  semicircle  approaches  zero  as  a— ^■oo.  Then 


J_ 

2k 


/ 


giK(z  -  z') 


dK 


± 


giK(z  -  z') 

^  -  k^ 

z 


dK. 


(89) 


By  the  residue  theorem,  the  integral  on  the  right  side  is  simply  2niE  (residues 
of  poles  inside  CJ,  which  in  the  limit  includes  all  poles  in  the  associated  half 
plane.  The  sign  on  the  right  side  of  equation  (89)  depends  on  which  contour 
is  used.  The  residue  theorem  requires  that  the  contour  be  traversed  in  a 
counterclockwise  direction.  If  Q  is  chosen  to  fall  in  the  lower  half  plane,  the 
interval  [-a,a\  is  traversed  from  a  to  -a,  so  the  negative  sign  must  be  used. 

To  avoid  poles  on  the  real  axis,  k,  is  taken  as  complex  with  imaginary  part 
6  >  0.  This  corresponds  to  absorption,  which  reduces  the  amplitude  of  a 
propagating  waveform  by  The  poles  of  the  integrand  occur  at  ±  k^  with 
one  in  the  upper  half  space  and  one  in  the  lower  half  space.  Thus,  only  one 
will  be  enclosed  by  Q.  For  the  case  z  -  z '  >  0,  the  contour  must  fall  in  the 
upper  half  plane,  and  the  resulting  residue  is 
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In  this  case,  the  contour  must  always  be  chosen  in  the  upper  half  plane  because 
the  argument  of  the  exponential  is  always  positive.  The  reflection  coefficient 
introduces  the  extra  pole  at  -p,  which  will  fall  in  the  upper  half  plane  because 
of  the  nature  of  Z„. 

O 

Using  equations  (92)  and  (93),  the  Green’s  function  may  be  written  as 
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(95) 
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the  integrations  over  height  z  and  vertical  wavenumber  k  necessary  for  forward 
and  inverse  Fourier  Transforms, 
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Exchanging  the  order  of  integration  gives  equation  (97)  for  the  homogenous 
propagation  term: 
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The  Fourier  Transforms  are  clearly  visible  in  the  first  two  lines.  The  first  line 
represents  the  direct  wave,  and  the  second  line  is  the  reflected  wave.  The  third 
line  represents  the  surface  or  creeping  wave  [12].  The  homogenous  phase 
multiplication  can  now  be  simplified  by  performing  the  contour  integral 
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(98) 
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The  residue  theorem  is  used  again,  and  the  poles  are  determined  by  recalling 
that  =  k/  -  k^.  The  denominator  is  then  -  (kg  -  t^),  with  zeros 

K  =  ±  W 

At  this  point,  it  is  necessary  to  determine  whether  the  poles  fall  in  the  spectrum 
of  Q.  For  the  moment,  assume  that  the  positive  pole  is  enclosed  by  the 
contour,  so 
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The  homogenous  propagation  term  then  becomes 
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(101) 


Compare  this  to  the  wide-angle  split-step  PE  and  recall  that  the  two  split-step 
approximations  <l>w  and  use  the  same  homogenous  propagation  term  A{ix), 
with  different  boundary  conditions.  As  mentioned  above,  the  wide-angle  PE 
implicitly  uses  a  pressure-release  boundary  condition.  In  terms  of  impedance, 
this  is  equivalent  to  Z^=0,  because  the  pressure  is  equal  to  zero.  In  this  case, 

R(k)  =  -1,  and  equation  (101)  agrees  exactly  with  the  homogenous  propagation 
part  of  the  wide-angle  split-step  PE.  Thus,  it  is  concluded  from  physical 
reasoning  that  equation  (100)  is  correct  and  equation  (101)  does  indeed 
reprs-  ent  the  homogenous  propagation  portion  of  the  GFPE. 
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Adding  the  nonhomogenous  phase  correction  term  yields  as  the  final  form  for 
the  GFPE 
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Note  the  two  differences  between  the  GFPE  and  wide-angle  PE.  One  is  the 
correction  factor  used  for  the  nonhomogenous  sound  speed  profile;  the  other  is 
the  reflection  coefficient.  The  wide-angle  PE  assumes  a  reflection  coefficient 
R  =  -1,  which  produces  no  amplitude  change  and  a  180°  phase  shift.  The 
GFPE  reflection  coefficient  produces  angular  dependant  amplitude  and  phase 
changes  and  also  introduces  the  surface  wave  term  that  results  from  the  zero 
ofR. 

In  practice,  the  first  two  lines  can  be  implemented  via  a  single  FFT/inverse  FFT 
pair.  Note  that  the  GFPE,  as  with  any  split-step  PE,  generates  a 
two-dimensional  (2-D)  pressure  field,  because  the  FFT  gives  a  profile  of 
pressure  versus  height  at  each  range  step.  It  is  generally  necessary  to 
interpolate  between  points  to  get  the  field  at  a  specified  height. 
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5.  Model  Comparisons 


The  FFP  [13]  is  another  model  in  wide  use  in  the  atmospheric  acoustics 
community.  It  is  a  full-wave  model,  based  on  an  integral  transform  of  the 
reduced  wave  equation.  The  FFP  transforms  with  respect  to  range  and 
horizontal  wavenumber,  thereby  preserving  the  boundary  condition  and 
height-dependant  sound-speed  in  the  transform  domain.  As  discussed  above, 
the  GFPE  is  an  integral  transform  in  height  and  vertical  wavenumber,  but 
incorporation  of  the  boundary  conditions  and  height-dependant  sound-speed  is 
not  straightforward. 

The  FFP  numerically  solves  for  a  height-dependant  Green’s  function  in  the 
horizontal-wavenumber  domain,  then  transforms  to  the  range  domain.  Because 
of  this,  the  sound-speed  profile  and  ground  impedance  must  be  range 
independent;  there  is  no  place  for  range  dependance  in  the 
horizontal-wavenumber  domain.  Thus,  it  exactly  incorporates  the 
height-dependant  sound-speed,  but  the  cost  is  range  independence. 

In  contrast,  the  Green’s  function  of  the  GFPE  represents  a  homogenous 
atmosphere  in  the  vertical-wavenumber  domain;  the  sound-speed  variations  are 
handled  by  an  approximation  in  the  spacial  domain.  The  advantage  is  the  sound 
speed  and  ground  impedance  can  vary  with  range. 

Figure  8  shows  the  full  2-D  output  of  the  GFPE  for  the  sound-speed  profile  of 
figure  7.  Note  that  this  profile  is  a  simulated  surface  duct,  not  a  measured 
profile.  Figure  9  gives  a  comparison  between  the  FFP  and  the  one-dimensional 
(1-D)  sound  pressure  level  versus  range  output  of  the  GFPE,  again  using  the 
profile  of  figure  7.  The  frequency  for  figures  8  and  9  is  20  Hz  and  the  source 
height  is  5  m.  The  blocky  appearance  of  the  2-D  20-Hz  output  is  due  to  the 
large  step  size  of  the  GFPE.  At  20  Hz,  with  a  reference  sound  speed  of  around 
340  m/s,  the  wavelength  is  about  17  m.  The  step  size  used  by  the  GFPE  was 
about  15  wavelengths,  which  equates  to  around  255  m.  At  low  frequencies, 
there  is  a  tradeoff  between  speed  and  resolution  in  the  GFPE. 

Figures  10  and  11  show  the  same  results  as  figures  8  and  9,  for  a  frequency  of 
100  Hz.  Figures  13  and  14  and  15  and  16  represent  the  GFPE  and  FFP  outputs 
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at  frequencies  of  20  and  100  Hz,  respectively,  for  the  upward-refracting  profile 
of  figure  12.  As  with  figure  7,  this  is  an  artificial  profile.  The  source  height 
is  again  5  m. 

The  profile  of  figure  17  represents  an  actual  surface  duct  calculated  from 
temperature  and  wind  measurements.  Figures  18  and  19  show  the  GFPE  output 
for  this  profile  at  20  and  100  Hz,  with  a  source  height  of  5  m.  Figure  20 
shows  another  profile  for  an  actual  surface  duct.  Figure  21  shows  the 
corresponding  GFPE  output  at  20  Hz.  Note  that  these  profiles  actually  contain 
multiple  ducts  near  the  surface,  giving  rise  to  a  fairly  complex  sound  field. 

Finally,  figure  22  shows  an  actual  upper-air  duct  profile  with  the  GFPE  output 
at  50  Hz.  In  figure  23,  the  source  height  is  500  m,  very  close  to  the  local 
minimum  in  the  profile.  This  tends  to  somewhat  focus  the  sound  energy.  Also 
note  the  shadow  zone  at  the  surface  beyond  8  km. 

As  mentioned  in  the  introduction,  a  parameter-sensitivity  study  was  conducted 
on  the  GFPE.  The  technical  report  for  this  study  contains  a  detailed  discussion 
of  the  input  parameters  and  a  much  more  extensive  presentation  of  the  GFPE 
output  and  comparisons  with  the  FFP. 
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Figure  8.  2-D  GFPE  output  for  the  sound  speed  profile  of  figure  7.  Frequency  is  20  Hz,  source  height  is 
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9.  Comparison  of  1-D  GFPE  and  FFP  outputs  for  the  sound  speed  profde  of  figure  7.  Frequency  is 
source  height  is  5  m. 
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Figure  10.  2-D  GFPE  output  for  the  sound  speed  profile  of  figure  7.  Frequency  is  100  Hz,  source  height  is 
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Figure  15.  2-D  GFPE  output  for  the  sound  speed  profile  of  figure  12.  Frequency  is  100  Hz,  source  height  is 
5  m. 
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e200910  Sound  Speed 


Measured  surface  duct  sound  speed  profile  associated  with  figures  18  and  19. 
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Figure  18.  2-D  GFPE  output  using  the  sound  speed  profile  of  figure  17.  Frequency  is  20  Hz,  source  height 
is  5  m. 
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Figure  19.  2-D  GFPE  output  for  the  sound  speed  profile  of  figure  17.  Frequency  is  100  Hz,  source  height  is 
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Measured  surface  duct  sound  speed  profile  associated  with  figure  21. 
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Figure  21.  2-D  GFPE  output  using  sound  speed  profile  of  flgure  20.  Frequency  is  20  Hz,  source  height  is 
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!  sound  speed  profile  of  figure  22.  Frequency  is  50  Hz,  source  height 
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Acronyms  and  Abbreviations 


FFP 

Fast  Field  Program 

FFT 

Fast  Fourier  Transform 

GFPE 

Green’s  function  parabolic  equation 

1-D 

one-dimensional 

PE 

parabolic  equation 

2-D 

two-dimensional 
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NAVAL  SURFACE  WEAPONS  CTR 
CODE  G63 

DAHLGREN  VA  22448-5000 

ARMY  OEC 

CSTE  EFS 

PARK  CENTER  IV 

4501  FORD  AVE 

ALEXANDRIA  VA  22302-1458 

ARMY  CORPS  OF  ENGRS 
ENGR  TOPOGRAPHICS  LAB 
ETL  GS  LB 

FT  BELVOIR  VA  22060 

TAC  DOWP 
LANGLEY  AFB 
VA  23665-5524 

ARMY  TOPO  ENGR  CTR 
CETEC  ZC  1 

FT  BELVOIR  VA  22060-5546 

LOGISTICS  CTR 
ATCL  CE 

FT  LEE  VA  23801-6000 

SCI  AND  TECHNOLOGY 
101  RESEARCH  DRIVE 
HAMPTON  VA  23666-1340 

ARMY  NUCLEAR  CML  AGCY 
MONA  ZB  BLDG  2073 
SPRINGFIELD  VA  22150-3198 

ARMY  FIELD  ARTLLRY  SCHOOL 
ATSF  F  FD 

FT  SILL  OK  73503-5600 


USATRADOC 
ATCD  FA 

FT  MONROE  VA  23651-5170 

ARMY  TRADOC  ANALYSIS  CTR 
ATRC  WSS  R 
WSMR  NM  88002-5502 

ARMY  RESEARCH  LABORATORY 
AMSRL  BE  M 
BATTLEFIELD  ENVIR  DIR 
WSMR  NM  88002-5501 

ARMY  RESEARCH  LABORATORY 
AMSRL  BE  A 

BATTLEFIELD  ENVIR  DIR 
WSMR  NM  88002-5501 

ARMY  RESEARCH  LABORATORY 
AMSRL  BE  W 
BATTLEFIELD  ENVIR  DIR 
WSMR  NM  88002-5501 

ARMY  RESEARCH  LABORATORY 
AMSRL  BE 
ATTN  MR  VEAZEY 
BATTLEFIELD  ENVIR  DIR 
WSMR  NM  88002-5501 

DEFNS  TECH  INFO  CTR 
CENTER  DTIC  BLS 
BLDG  5  CAMERON  STATION 
ALEXANDRIA 
VA  22304-6145 

ARMY  MISSILE  CMND 
AMSMI 

REDSTONE  ARSENAL 
AL  35898-5243 

ARMY  DUGWAY  PROVING  GRD 
STEDP  3 

DUGWAY  UT  84022-5000 

USATRADOC 
ATCD  FA 

FT  MONROE  VA  23651-5170 


ARMY  FIELD  ARTLRY  SCHOOL 
ATSF 

FT  SILL  OK  73503-5600 

WSMR  TECH  LIBRARY  BR 
STEWS  IM  IT 
WSMR  NM  88001 


Record  Copy 


